/*  sdfasins.cpp  by K.Tsuru  */
// function ID 3204   DRADIX
/*****************************************************************
SDouble class
inverse trigonometric function arcsin x by series for |x| < 1.0
If fe <= 0, its value is taken as the value of fixed exponent.
Default argument is declared as fe = INT_MAX > DRADIX_EXP_MAX.
If |x| << |a| and you would like to get an addition
    y = a + AsinSeries(x),
it may be better to modify into
    y = a + AsinSeries(x, a.RdxExp()).
*****************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static const char* func = "AsinSeries";
SDouble AsinSeries(const SDouble& x, int fe){
    if(x.NetRdxExp() > 0) x.SetError(x.DOMAIN_ERR, func, 3204); // |x| > 1.0
    SDouble sum(x), xsq, delta;
    ulong n2 = 4, den, num, mt = sum.SlOpMaxValue();

    fe = (fe > DRADIX_EXP_MAX) ? x.RdxExp() : fe;
    sum.FixedPoint(fe);
    xsq = x*x;
    delta = DsDiv(x*xsq, 6);
    sum += delta;
    do{
        num = (n2-1)*(n2-1);
        den = n2*(n2+1);
        if(den >= mt)   break;
        delta = DsMult(delta, num);
        delta = DsDiv(delta*xsq, den);
        sum += delta;
        n2 += 2;
    } while(delta.Sign(3204));

    while(delta.Sign(-3204)){
        if(n2 + 1 >= mt){
            sum.SetError(sum.NOT_CONVERGE, func, -3204);
            break;
        }
        delta = DsMult(delta*xsq, (n2-1));
        delta = DsMult(delta, (n2-1));
        delta = DsDiv(delta ,n2);
        delta = DsDiv(delta ,n2+1);
        sum += delta;
        n2 += 2;
    }
    sum.PointFree();
    sum.Reform(3204);
    sum.upToTerm = n2/2;

    return sum;
}